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Abstract 

Structural reliability methods aim at computing the probability of failure of systems with 
respect to some prescribed performance functions. In modern engineering such functions 

usually resort to running an expensive-to-evaluate computational model (e.g. a finite el- 
ement model). In this respect simulation methods, which may require 10''"^ runs cannot 
be used directly. Surrogate models such as quadratic response surfaces, polynomial chaos 
expansions or kriging (which are built from a limited number of runs of the original model) 
are then introduced as a substitute of the original model to cope with the computational 
cost. In practice it is almost impossible to quantify the error made by this substitution 
though. In this paper we propose to use a kriging surrogate of the performance function as 
a means to build a quasi-optimal importance sampling density. The probability of failure is 
eventually obtained as the product of an augmented probability computed by substituting 
the meta-model for the original performance function and a correction term which ensures 
that there is no bias in the estimation even if the meta-model is not fully accurate. The ap- 
proach is applied to analytical and finite element reliability problems and proves efficient 
up to 100 random variables. 

Keywords: reliability analysis, importance sampling, metamodeling error, kriging, random 
fields, active learning, rare events 



1, Introduction 

Reliability analysis consists in the assessment of the level of safety of a system. Given 
a probabilistic model (an n-dimensional random vector X with probability density function 
(PDF) fx) and a performance model (a function g), it makes use of mathematical tech- 
niques in order to estimate the safety level of the system in the form of a failure probability. 
A basic reference approach is the Monte Carlo simulation technique that resorts to numeri- 
cal simulation of the performance function through the probabilistic model. 

Failure is usually defined as the event F = {g (X) < O}, so that the failure probabihty is 
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defined as follows: 



p^=P(F) = P({g(X)<0}) 

fxMdx 

% = {x6R" :g(x)<0} 

Introducing the failure indicator function lg<o being equal to one if g (x) < and zero oth- 
erwise, the failure probability turns out to be the mathematical expectation of this indicator 
function with respect to the joint probability density function fx of the random vector X. 
This convenient definition allows one to derive the Monte Carlo estimator which reads: 

k=i 

where ^x^^\ ... ,x^'^^, AT e N* is a set of samples from the random vector X. Accord- 
ing to the central limit theorem, this estimator is asymptotically unbiased and normally 
distributed with variance: 

In practice, this variance is compared to the unbiased estimate of the failure probability in 
order to decide whether the simulation is accurate enough or not. To do so, one usually re- 
sorts to the coefficient of variation defined as 5^^^^^ = crp^^^/P/MC- ^ being fixed, this coeffi- 
cient of variation dramatically increases when the failure event becomes rare (pf 0). This 
makes the Monte Carlo estimation technique intractable for real world engineering prob- 
lems for which the performance function involves the output of an expensive-to-evaluate 
black-box function - e.g. a finite element code. Note that this remark is also true for too 
frequent events (p^ — > 1) as one should compute the coefficient of variation of 1 — P/mc 
that exhibits the same property. 

In order to reduce the number of simulation runs, different alternatives to the brute- 
force Monte Carlo method have been proposed and might be classified as follows. 

One first approach consists in replacing the original experiment by a surrogate which is 
much faster to evaluate. Various surrogates have been used amongst which are: quadratic 
response surfaces (Bucher and Bourgund, 1990; Kim and Na, 1997; Das and Zheng, 2000) 
and the more recent metamodels such as support vector machines (Hurtado, 2004; Deheeger 
and Lemaire, 2007), neural networks (Papadrakakis and Lagaros, 2002) and kriging (Kay- 
maz, 2005; Bichon et al., 2008). Nevertheless, it is often difficult or even impossible to 
quantify the error made by such a substitution. 

The well-known first- and second-order reliability methods, based on Taylor expan- 
sions somewhat differ from these surrogate-based approaches because of their mathemat- 
ical background (Breitung, 1984). But in practice they feature the same limitation: the 
assumptions (mostly the unicity of the so-called most probable failure point) they are built 
on may not hold and it is difficult to validate them. 
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Bucher (2009); Naess et al. (2009); Nishijima et al. (2010) investigate the use of the 
extreme value theory in order to infer the tail of the cumulative distribution function (CDF) of 
the random variable g(X) from a set of samples. This set is obtained by crude Monte Carlo 
simulation of the input random vector X propagated through the performance function g. 
Despite this method has a lot of attractive features (it is independent of both the dimension 
n and the shape of g), it may require a large number of samples to make the fitted extreme 
value distribution accurate far in the tail - i.e. for low failure probabilities. 

From another point of view, variance reduction techniques have been proposed in or- 
der to make Monte Carlo simulation more efficient. Importance sampling (Rubinstein and 
Kroese, 2008) aims at concentrating the Monte Carlo samples in the vicinity of the limit- 
state surface, e.g. around the most probable failure point (also known as design point) 
obtained by a preliminary FORM analysis (Melchers, 1989). Subset simulation (Au and 
Beck, 2001; Ching et al., 2005a,b; Santoso et al., 2010) computes the failure probability 
as a product of conditional probabilities, each of them being estimated by Markov Chain 
Monte Carlo simulation. All these approaches reveal robust, although they are often too 
much computationally demanding to be implemented in industrial cases. 

As a summary the current practice for evaluating probabilities of failure in case of com- 
putationally demanding performance functions relies on the substitution of the limit-state 
function by a metamodel for which no general error estimation is usually available. 

In this paper, a new hybrid approach combining importance sampling and an adaptive 
metamodeling technique is proposed. First, a kriging surrogate of the limit-state function 
is built and adaptively refined. Then, the probabilistic prediction provided by the kriging 
surrogate is used to build up a quasi-optimal importance sampling density. As a result 
the probability of failure is computed as the product of two terms, one which is estimated 
by sampling the surrogate limit-state function, and a correction factor computed from the 
original limit-state function. 

The paper is organized as follows. Section 2 recalls the basic of Gaussian process 
(kriging) metamodeling and introduces the probabilistic classification function. Section 
3 presents the construction of a quasi-optimal importance sampling density derived from 
this function and the associated estimator of the failure probability. Section 4 introduces an 
adaptive refinement technique of this importance sampling density function so as to make 
the algorithm as parsimonious as possible. Section 5 and 6 eventually describe practical 
implementation details and application examples. 

2. Probabilistic classification using margin metamodels 

A metamodel means to a model what the model itself means to the real-world. Loosely 
speaking, it is the model of the model. As opposed to the model, its construction does not 
rely on any physical assumption about the phenomenon of interest but rather on statistical 
considerations about the coherence of some scattered observations that result from a set 
of experiments. This set is usually referred to as a design of experiments (DOE): 3^ — 
{xi, . . . ,x^}. It should be carefully selected in order to retrieve the largest amount of 
statistical information about the underlying functional relationship over the input space 
S!^. In the sequel one attempts to build a metamodel for the failure indicator function lg<o- 
In the statistical learning theory (Vapnik, 1995) this is referred to as a classification problem. 
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In this paper, we define a margin metamodel as a metamodel that is able to provide 
a probabilistic prediction of the response quantity of interest whose spread (i.e. variance) 
depends on the amount of information brought by the DOE. Note that this spread is thus 
reducible by bringing more observations into the DOE. In other words, this is an epistemic 
(reducible) source of uncertainty that shall not be confused with the possible uncertainty 
in X that will be modelled by a PDF fx later on for structural reliability purposes. To 
the authors' knowledge, there exist only two families of such margin metamodels: the 
probabilistic support vector machines P-SVM by Piatt (1999) and Gaussian process- (or 
kriging-) based classification. The present paper makes use of the kriging metamodel, 
but the overall concept could easily be extended to P-SVM. The theoretical aspects of the 
kriging prediction are briefly introduced in the following subsection before it is applied to 
the classification problem of interest. 

2.1. Gaussian-process based prediction 

The Gaussian process (also known as kriging) metamodeling theory is detailed in Sacks 
et al. (1989); Welch et al. (1992); Santner et al. (2003). In essence, it assumes that the 
performance function g is a sample path of an underlying Gaussian stochastic process G 
that may be cast as follows: 

GM^fixf P+Zix) (4) 

where: 

• f (x)^ p denotes the mean of the GP which corresponds to a classical linear regres- 
sion model on a given functional basis {f^, i = 1, . . . , p} G ^2 i^x> ^] ^ 

• Z (x) denotes a zero-mean stationary Gaussian process with a constant variance cjg. 
It is fully defined by its autocovariance function which reads: 

CaGix,x')=cjlR{x-x\e) (5) 

where £ is a vector of parameters defining R. 

The most widely used class of autocorrelation functions is the anisotropic squared exponen- 
tial model: 



n 



R(.-x',0=exp(-l](^i^^j 1 (6) 

The best linear unbiased estimation (BLUE) of G at point x is shown (Santner et al., 
2003) to be the following Gaussian random variate: 



2^ 



G(x) = Gix) |{G(Xi) = g(Xi), . . .,GixJ = gixj} 
~^(/ig(x), (jg(A:)) 



(7) 



with moments: 



^JiG(x)^f(xyp + r(xfR-'[y-¥p) (8) 
cr|(x) = cr2 (l-r(x)"^R-V(x) + u(x)^(F^R-^F)-^u(x)) (9) 
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where y = (gix^), ... , g{x^) is the vector of observations, R is their correlation matrix, 
r(x) is the vector of cross-correlations between the observations and the prediction, and F 
is the so-called regression matrix. The terms of R, r{x) and F are respectively defined as 
follows: 

ri{x)=R{x-x„t),i^l,...,m (10) 
R,j^R(x,-Xj,l), i, j = l, ...,m (11) 
^0 =/j (^i)' i = 1, ■■■ ,m, j = 1, ... ,p (12) 

Finally, the generalized least-squares solution P and the vector u{x^ respectively read: 

^ = (F"^R-^F)-1f"^R-V (13) 
u(x) = F"^R-^r(x)-/(x) (14) 

At this stage one can easily prove that jUg (xj = g (x;) and erg (xj = for i = 1, . . . , m, 
thus meaning the kriging surrogate interpolates the observations without any residual epis- 
temic uncertainty. 

Given a choice for the regression and correlation models, the optimal set of parameters 
, C and CTg* can then be inferred using the maximum likelihood principle applied to the 
single sparse observation of the GP sample path grouped into the vector y. This inference 
problem turns into an optimization problem that can be solved analytically for both fi* 
and CTg* assuming is known. Thus the problem is solved in two steps: the maximum 
likelihood estimation of is first solved by a global optimization algorithm which in turns 
allows one to evaluate the optimal ^* and cr^*. Implementation details can be found in 
Welch et al. (1992); Lophaven et al. (2002). 

2.2. Probabilistic classification function 

As seen from Eq. (7), the kriging surrogate provides both an approximation of the limit- 
state function g(x) which is denoted by /ig(x) and an epistemic prediction uncertainty 
which is characterized by the kriging variance a^{x). 

Authors usually make use only of the surrogate g{x) = /ig(x) for reliability analysis - 
see Kaymaz (2005); Bichon et al. (2008). In this paper it is proposed to use the complete 
probabilistic prediction - see Picheny (2009) for a similar idea. 

Let us introduce for this purpose the following probabilistic classification function: 

7i(x) = ^ [g(x) < O] (15) 

In this expression, the probability measure ^ [•] refers to the Gaussian nature of the krig- 
ing predictor G{x) ~ ^(jLtg(x), crg(x)) and shall not be confused with the probability 
measure P(») associated with the random vector X - see Eq. (1). 

Thanks to the Gaussian nature of the kriging predictor, the probabilistic classification 
function rewrites: 

/^0-Ltg(x)^ 

7r(x) = * if;c^^ (16) 
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and still holds for the points in the experimental design for which the prediction variance 
equals zero by switching to the limit: 



/ 1 if XG^,gix)<0 

It shall be again emphasized that 7i(x) is not the sought failure probability. It may be 
interpreted as the probability that the predictor G(x) (for some prescribed deterministic x) 
is negative with respect to the epistemic uncertainty. 

Figure 1 illustrates the concepts introduced in this section on a basic structural reli- 
ability example from Der Kiureghian and Dakessian (1998). This example involves two 
independent standard Gaussian random variates X-^ and X2, and the performance function 
reads: 

g(xi, X2) = b -X2 -?c(xi -e)^ (18) 

where b — 5, k — 0.5 and e = 0.1. In Figure 1, the original limit-state function gix) — 
is represented by the black solid line. The red and the blue "— " represent the initial 
DOE from which the kriging metamodel is built, the color being related to the sign of the 
performance function g. Once kriging has been applied, the mean prediction of the limit- 
state surface, which is defined as {x e M" : jJigix) — O}, is represented by the dashed black 
line. It can be seen that the metamodel is not fully accurate since the green triangle x° 
(among others) is misclassified. Indeed, x° is in the safe domain according to the real 
performance function g, and in the failure domain according to the mean kriging predictor 

Note that the probabilistic classification function allows a smoother decision: x° belongs 
to the failure domain with a 60% probability wrt. the epistemic uncertainty in the random 
prediction G(x°) ~ ^(/>ig(x°), crg(x°)). Note also that the red and the blue "-" in the 
DOE are always surely classified due to the interpolating property of the kriging metamodel 
= Ptg(^i), i = 1, ... ,m). 

Figure 2 is the one-dimensional illustration of the two classification strategies at point 
x° (green triangle): 



according to the deterministic decision function (the Heaviside function centered in 
zero), the real performance function g is positive whereas the mean kriging prediction 
jJLQ is negative; 

according to the probabihstic classification function (the smoother Gaussian cumu- 
lative distribution function centered in zero and whose spread is characterized by 
crg(x°)), the epistemic prediction G is negative with a 60% probability (7i(x°) = 0.6 
as defined from Eq. (16)). 



3. Metamodel-based importance sampling 

Picheny (2009) proposes to use the probabilistic classification function n (see Eq. (16)) 
as a surrogate for the original indicator function l(,<o, so that the failure probability is 
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Figure 1 : Comparison of the two classification strategies on a two-dimensional example from Der Kiureghian 
and Dakessian (1998) - Limit-state functions and kriging probability values. 




Figure 2: Comparison of the two classification strategies on a two-dimensional example from Der Kiureghian 
and Dakessian (1998) - Classification functions at point 
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rewritten from its definition in Eq. (1) as follows: 



Pfs = 



7i(x)fY(x)dx 

jjxy ^^^^ 



It is argued here that this latter quantity does not equal the failure probability of interest 
because it sums the aleatory uncertainty in the random vector X and the epistemic uncer- 
tainty in the prediction G. This is the reason why p^^ will be referred to as the augmented 
failure probability in the sequel. As a matter of fact, even if the epistemic uncertainty in 
the prediction can be reduced (e.g. by enriching the DOE as proposed in Section 4), it is 
impossible to quantify the contribution of each source of uncertainty in p^^. 

This remark motivates the approach introduced in this section where the probabilistic 
classification function is used in conjunction with the importance sampling technique in 
order to build a new estimator of the failure probability. 

3. 1 . Importance sampling 

According to Rubinstein and Kroese (2008), importance sampling (IS) is one of the 
most efficient variance reduction techniques. This technique consists in computing the 
mathematical expectation of the failure indicator function according to a biased PDF which 
favors the failure event of interest. This PDF is called the instrumental density. Let h denote 
such a probability density such that h dominates lg<o/xj meaning that: 

Vxe®^, h(x) = 0^1g<o(A:)/;f(A:) = (20) 

Given this instrumental density, the definition of the failure probability of Eq. (1) may 
be rewritten as follows: 



Pf 



ig<o(^)i:r^K:>c)dx 

h(x) 

(21) 

fxixy 



In this expression, the expectation is now computed with respect to the instrumental density 
h. 

The latter definition of the failure probability easily leads to the definition of the impor- 
tance sampling estimator: 



1 ^ fxix^^^) 



(22) 



where {x^^^ . . . ,x^~)}, AT e N* is a set of samples drawn from the instrumental density h. 
According to the central limit theorem, this estimation is unbiased and its quality may be 
measured by means of its variance of estimation which reads: 

1 f 1 ^ ffx^^^l^ A 
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Rubinstein and Kroese (2008) show that this variance is zero (optimaUty of the IS esti- 
mator) when the instrumental PDF is chosen as: 



h*(x) = 



l,<oix)fix) 



j lg<o(x)/(x)dA: 
lg<o(x)/(x) 

Pf 



(24) 



However this instrumental PDF is not implementable in practice because it involves the 
sought failure probability in its denominator. The art of importance sampling consists in 
building an instrumental density which is quasi-optimal. 

3.2. A metamodel-based approximation of the optimal instrumental PDF 

Different strategies have been proposed in order to build quasi-optimal instrumental 
PDF suited for specific estimation problems. For instance, Melchers (1989) uses a standard 
normal PDF centered onto the design point obtained by FORM in the space of the indepen- 
dent standardized random variables. This approach might be inaccurate if the design point 
is not unique though. Cannamela et al. (2008) use a kriging prediction of the performance 
function g in order to build an instrumental PDF suited for the estimation of extreme quan- 
tiles of the random variate g{X). However this approach is not applied to the estimation of 
failure probabilities. 

In this paper it is proposed to use the probabilistic classification function in Eq. (16) as 
a surrogate for the real indicator function in the optimal instrumental PDF in Eq. (24). The 
proposed quasi-optimal PDF thus reads as follows: 



where ^ is the augmented failure probability which has been already defined in Eq. (19). 
For the sake of illustration, this quasi-optimal instrumental PDF is compared to the optimal 
(although impractical) one in Figure 3 using the example of Section 2.2. 

3.3. The metamodel-based importance sampling estimator 

Choosing the proposed quasi-optimal instrumental PDF (i.e. substituting h* for h in 
Eq. (21)) leads to the following new expression of the failure probability: 



h*(x) = 



7I(X)/(X) 




(25) 



<x)fix) 

Pfs 




r 



Pf = lg<o(^) 



h*ix)dx 



(26) 



h*{x)dx 



(27) 



corr 



(28) 
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(a) The optimal instrumental PDF h*. 



(b) The proposed quasi-optimal instrumental PDF h*. 



Figure 3: Comparison of the instrumental PDFs on the two-dimensional example from Der Kiureghian and 
Dakessian (1998). 



where we have introduced: 



^corr — '^h* 



lg<o(X) 



(29) 



This means that the failure probabihty is equal to the product of the augmented failure 
probability p^g and a correction factor a^^^^. 

This correction factor is defined as the expected ratio between the real indicator function 
lg<o and the probabilistic classification function n. Thus, if the kriging prediction is fully 
accurate, the correction factor is equal to one and the failure probability is identical to 
the augmented failure probability (optimality of the proposed estimator). On the other 
hand, in the more general case where the kriging prediction is not fully accurate (since it 
is obtained from a DOE 3^ of finite size m), the correction factor modifies the augmented 
failure probability accounting for the epistemic uncertainty in the prediction. 

The two terms in the latter definition of the failure probability may now be estimated 
using Monte Carlo simulation: 



i=l 



N. 



(30) 
(31) 



where the first N^-sample set is generated from the original PDF fx, and the second N^on' 
sample set is generated from the quasi-optimal instrumental PDF h*. According to the 
central limit theorem, these two estimates are unbiased and normally distributed. Their 
respective variance denoted by cr^ and 0"^^^.^. may be easily derived as in Eq. (23). 
Finally, the proposed estimator of the failure probability simply reads as follows: 



Pf metalS ~ Pf e ^corr 



(32) 
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It is important to note that both terms in Eq. (32) are independent, since they rely upon 
samphng according to two independent PDFs. Based on this latter remark, it is shown in 
Appendix A that for reasonably small values of the coefficients of variation 5^^^ and 5^^^^ of 
the two estimators p^r^ and a^^^ (say 1 - 10%), the coefficient of variation of the proposed 
estimator approximately reads: 

5g ^ J 5^ +51 (33) 

3.4. Sampling from the quasi-optimal instrumental density 

The instrumental density h* in Eq. (25) is not straightforward to sample from. First one 
observes that it is defined up to an unknown but finite normalizing constant: 

h*(x)oc7i(jc)/(x) (34) 

Indeed the exact value of p^^ is not known (it will only be estimated by Monte Carlo simu- 
lation). It is clear that the usual inverse transform simulation techniques are not applicable 
for this non-parameteric multidimensional pseudo-PDF. This is the reason why one resorts 
to Markov chain Monte Carlo (MCMC) simulation techniques (Robert and Casella, 2004). 

These techniques consists in generating a first-order Markov chain 2^ — jz*^'-*, i = 1, . . . ,N 
whose asymptotic distribution is shown to be the target distribution (see Robert and Casella, 
2004, pp. 206 - 207). The present implementation makes use of the slice sampling tech- 
nique (Neal, 2003) - see Appendix B. 

However, Robert and Casella (2004) also point out that the Markov chain samples are 
not independent and identically distributed: there is some correlation between the suc- 
cessive states of the chain. Thus, in order to ensure that the samples used in Eq. (30) 
are independent and identically distributed, we use the so-called thining procedure which 
consists in keeping only one sample every t states of the chain (say t — 10). 



4. Adaptive refinement of the probabilistic classification function 

The significance of the variance reduction introduced by the proposed metamodel-based 
importance sampling technique mostly relies on the optimality of the proposed instrumental 
PDF h*. Thus, it is proposed here to adaptively refine the probabilistic classification function 
so that the quasi-optimal instrumental PDF h* converges towards its optimal counterpart 
h*. 

4. 1 . Refinement strategy 
4.1.1. A short state-of-the-art 

There exists a relatively large literature about the adaptive refinement of a kriging pre- 
diction for accurate classification (or level-set approximation): see e.g. Oakley (2004); Lee 
and Jung (2008); Bichon et al. (2008); Vazquez and Beet (2009); Picheny et al. (2010); 
Echard et al. (2011). They all rely on the definition of a so-called in-fill criterion which is 
maximum or minimum in the region of interest: the region where the sign of the predicted 
performance function is the most uncertain. The interested reader may find the expressions 
for all these criteria together with a discussed comparison on two analytical examples in a 
recent paper by Beet et al. (2011). 
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For instance, Echard et al. (2011) propose the so-called -criterion which is defined as 
follows: 



(35) 



crg(x) 

This function is minimum in the vicinity of the limit-state surface of the mean prediction 
{x e E" : fiQix) — O}. This vicinity corresponds to the grayest part of the shade area in 
Figure 1. 

Then the problem of finding the best improvement point is usually reduced to a global 
optimization problem. Using the '^-criterion from Echard et al. (2011), the problem reads 
as follows: 

X = arg min '^{x) (36) 



'X 



and is usually solved by means of a numerical global optimization algorithm. 

Additionally it is argued that for the present reliability estimation problem, one should 
focus on the most "probabilistically significant" region - such as it is done in FORM, because 
this is the region where the misclassified points have the most significant impact on the 
optimality of the proposed estimator. In Echard et al. (2011), such a trade-off is achieved by 
picking the best improvement point in a discrete population X that is distributed according 
to fx- The global optimization is thus reduced to a discrete optimization problem: 

x^+i = argmin'^(x) (37) 

4. 1 .2. The proposed strategy 

As an introduction, it is argued that the available approximation of the optimal instru- 
mental PDF after m observations: 

/?(x)oc7i(x)/x(x) (38) 

is an interesting trade-off criterion because it is maximum (i) where the sign of the pre- 
dicted performance function is supposed to be negative, and (ii) where the probability 
density function fx is maximum. 

Then, it is proposed to use a sampling alternative to the optimization step - as in 
Dubourg et al. (2011). This is achieved by considering the in-fill criteria as PDFs defined 
up to an unknown but finite normalizing constant. It proceeds as follows: 

1. Sample a large population (say N — samples) from the weighted in-fill criterion 
(here h*) using a MCMC simulation technique such as slice sampling (Neal, 2003). 

2. Reduce this population to its K clusters' center (K being given, see Section 6) using 
the iiC-means clustering algorithm (MacQueen, 1967). 

3. Evaluate the performance function on the K clusters' center. 

4. Enrich the former experimental design with these K new observations. 

5. Update the kriging prediction and loop back to step 1 until some target accuracy is 
achieved - see Section 4.2. 

The proposed samphng-based refinement strategy allows one to refine the prediction 
from a batch of optimal points instead of a single best point. It thus solves the problem 
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of locally optimal points. Indeed, the in-fiU criteria proposed in the literature commonly 
features several optima thus meaning that there does not exist a single best point. In the 
proposed strategy all the maxima are interpreted as modes of the improper PDF and results 
in local concentrations of points in the first step. The K-means clustering algorithm used 
in the second step reduces the population generated in the first step to the most signifi- 
cant modes. Note also that this approach allows one to perform several computations of 
the performance model in a distributed manner (i.e. on a high performance computing 
platform) . 

4.2. Stopping criterion 

Finally a metric is required to check the optimality of the probabilistic classification 
function. This metric may then be used as a stopping criterion for the previously detailed 
refinement strategy thus making it adaptive. The metric is based on a cross-validation 
procedure. 

Cross-validation techniques are classically used for model selection (see e.g. Stone, 
1974). Basically, the design of experiments = {x^,i — 1, . . . , m} is split into a learning 
subset 3^^ and a validation subset 3^y such that n = and X^U 3^y ^ 3^. The 
model is then built using the learning subset S^i^ (hence denoted by and validated by 
comparing the predicted values g^^-^ix) and the real values g(x) onto the validation subset 
X G 3Cy. The leave-one-out technique is a special case where the learning subset is defined 
diS = X \ x^. In a regression context, Allen (1971) propose to use a leave-one-out 
estimate of the mean squared error referred to as the predicted residual sum of squares: 

1 

PRESS = -Xi (^^\x.(^.) - gi^i^y (39) 

i=l 

The metric proposed in this paper is a leave-one-out estimate of the correction factor in 
Eq. (29): 

where G^\^. is the i-th leave-one-out kriging prediction of the performance function built 
from the DOE without the i-th sample X;. This quantity should be estimated using a 
minimal number of samples (say m = 30) so that it is sufficiently accurate. 

The reason for introducing this leave-one-out estimate of a^^rr is the following. In the 
early steps of the refinement procedure, the kriging predictor is not accurate enough to 
evaluate the failure probability by means of Eq. (32) efficiently. It would thus be ineffective 
to waste costly evaluations of the true limit-state function to compute the true correction 
factor in Eq. (31). On the contrary, the leave-one-out estimate ScorrLoo makes only use of 
the available observations in the experimental design, so that it is fast to evaluate. As dis- 
cussed earlier in Section 3.3, a sufficient accuracy is reached when the estimate in Eq. (40) 
gets close to one because it means that the probabilistic classification function converges to- 
wards its deterministic counterpart. And consequently the instrumental PDF h* converges 
towards it optimal counterpart h*. It is thus proposed to use acorrLoo an indication to 
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stop the refinement strategy. Precisely, the iterative enrichment of X is stopped if its size is 
greater than 30 and Scoj-rLoo is in the order of magnitude of 1. In the case of high dimen- 
sional and/or highly nonlinear problems, the size of the DOE is limited to ra^^ = 1000. 
One may also measure the improvement brought by the latest refinement iteration in terms 
of variation of S^q^loo- 

Remark regarding numerical implementation. It may happen that a leave-one-out estimate 
of the kriging variance cr5. (xj gets close or even equal to zero. This is problematic 

because the ratio in Eq. (40) might not be defined in such cases. Indeed, if the mean 
prediction i^tg^^^ (^J is positive, the probabilistic classification function in the denominator 
equals zero - and it may cause an exception. It is argued that this variance should never be 
exactly zero (the kriging variance equals zero only at the samples in the DOE, here SCXx^), 
so that it is proposed to bound the probabilistic classification function above a reasonably 
low value (say the machine precision, % ^ 10"^^). 

5. Algorithmic implementation 

The purpose of this section is to summarize the proposed metamodel-based importance 
sampling algorithm. The flowchart of the algorithm is given in Figure 5. It is essentially 
divided in three independent steps, two of which could potentially be run in parallel. 

First, the algorithm only requires the choice of a maximum number of performance 
function evaluation N^^-^, a targeted coefficient of variation 5^^^^^^ for the final estimate 
P/metais ^^e number K of points that should be added at each refinement iteration. 

Then, the first step of the algorithm is to build a reasonably accurate approximation 
of the probabilistic classification function that will bring a significant variance reduction 
for the final importance sampling estimate. Note that this step is optional in the case the 
analysis is resumed from an existing DOE. 

The two other steps are independent and might be run in parallel. They consist in us- 
ing Monte Carlo simulation for the estimation of the two terms a^^^ and ^ defining the 
proposed estimator P/metais ~ see Section 3. In the current implementation, we targeted the 
same coefficient of variation 5^ — S^^^ = Sjarget/v^ for the two estimators. It could also 
be possible to implement these two steps in two independent threads that would commu- 
nicate between each other and decide whether to stop the simulation or not based on the 
estimate of the final coefficient of variation. It should also be pointed out that it is easier to 
reduce the estimation variance on the augmented failure probability than the one one the 
correction factor because the former only depends on the kriging surrogate. 

6. Application examples 

In this section, the proposed strategy is applied to two structural reliability problems. 
The first one involves a simple analytical limit-state function and may thus be used to val- 
idate the implementation of the proposed strategy whereas the second example involves 
a more sophisticated nonlinear finite element model coupled with a high dimensional 
stochastic model so as to prove the applicability of the strategy to industrial problems. 
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Figure 4: Flowchart of the proposed metamodel-based importance sampling algorithm 
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6.1. Analytical limit-state function - Influence of the dimension 

This basic structural reliability example is taken from Rackwitz (2001). It involves n 
independent lognormal random variates with mean value ji — 1 and standard deviation 
cr = 0.2. The performance function reads as follows: 

n 

= (n + ac^^/r^) -^X; (41) 

1=1 

where a is set equal to 3 for the present application. The dimension of the problem is 
successively set equal to n = {2, 50, 100} to assess the influence of the dimension on the 
proposed estimator. The results are given in Table 1. 



n 


2 


50 


100 


Crude Monte Carlo (ref.) 


Pf MC 


4.78 X 10"^ 


1.91 X 10-3 


1.73 X 10-3 




< 2% 


< 2% 


< 2% 


N 


522,000 


1,100,000 


1,450,000 


Metamodel-based importance sampling 


■'^DOE 


6x2 


6 X 50 


6 X 100 




5.03 X 10-3 


1.95 X 10-3 


1.83 X 10-3 




< 1.41% 


< 1.41% 


< 1.41% 


N 

^ ^ corr 


100 


1,500 


2,100 


'^corr 


1 


0.99 


0.93 


^corr 


0% 


< 1.41% 


< 1.41% 


^DOE + ^fcorr 


112 


1,800 


2,700 


Pf metalS 


5.03 X 10-3 


1.93 X 10-3 


1.70 X 10-3 


^metalS 


< 1.41% 


< 2% 


< 2% 



Table 1: Comparative results for the performance function g(x) = [n + auy/n) — from Rackwitz 

z"=l 

(2001). 

The first block provides the results from a crude Monte Carlo simulation (reference) 
whereas the second block provides the results from the proposed importance sampling 
scheme. For both methods a 2% coefficient of variation on the failure probability estimates 
is targeted. First, it can be observed that the proposed estimator is in good agreement 
with the reference Monte Carlo estimates. The estimates of both the augmented failure 
probability g and the correction factor a^^^j. are also provided. Note that the correction 
factor increases with the dimension. In low dimension (n = 2) the kriging surrogate almost 
exactly equals the performance function (at least in the vicinity of its limit-state) hence 
^corr — 1 ^i^d P/metais ~PfE (^^^ misclassificatiou) . In larger dimensions the surrogate loses 
accuracy, and the correction factor gets more and more influent. 

The size of the DOE from which the kriging prediction is built is given as the product 
between the number of refinement iterations and the number K of clusters' center added 
per iteration. K is chosen equal to the number of input random variables (n = {2, 50, 100}) 
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Figure 5: Illustration of the Scordelis-Lo shell roof example. 



SO that the cost for the definition of the instrumental PDF is equivalent to the cost induced 
by a gradient-based design point search algorithm (when the gradient is approximated 
using a finite forward dijferentiation technique). In all cases 6 iterations are required (see 
Table 1). Note that the total number of calls to the real performance function g for the 
proposed importance sampling scheme (Nqoe + ^corr) is much less than the number of calls 
induced by Monte Carlo simulation (N) for the same targeted coefficient of variation. 

6.2. Nonlinear stability analysis of an imperfect shell roof 

The mechanical model for this example is inspired by Scordelis and Lo (1961). It con- 
cerns the buckling analysis of a cylindrical shell roof whose dimensions are given in Figure 
5. The longitudinal edges of the roof are free while its circumferential edges are supported 
by rigid diaphragms (radial displacement fixed to zero). Its constitutive material is as- 
sumed to have a nonlinear elastic Ramberg-Osgood material behavior. It is subjected to a 
constant surface load q and the structure is considered to fail if its critical buckling load q^r 
is less than a prescribed service load of magnitude — 0.18 MPa, so that the associated 
performance function may be defined as follows: 

g(^) = qJ?)-qo (42) 

where ^ denotes the outcome of the random vector S introduced in the sequel. 

The critical buckling load q^r is determined by means of the asymptotic numerical method 
(Cochelin, 1994) coupled with a 30 x 30 8-node Biichter-Ramm shell finite element mesh 
(Biichter et al., 1994) using the EVE finite element code. This academic software was 
initially developed by Cochelin (1994) and further developed by Noirfalise et al. (2008). 
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Method Subset (ref.) Multi-FORM Meta-IS 



DOE size 



6 X 93 



MPFP search 



10, 000 



Simulations 



20,000 
1.27 X lO""* 



9,464 
1.32 X 10-"^ 



Pf estimates 



1.22 X lO-'^ 



C.o.V 



12.36% 



13.75% 



Table 2: Comparative results for the shell roof example adapted from Dubourg et al. (2009). 

The Stochastic model, inspired from Dubourg et al. (2009), involves four independent 
random fields defined over the roof surface. They describe respectively: the initial shape 
imperfection the material Young's modulus E, the material yield stress cr^, and the shell 
thickness h. The random shape imperfection is modeled as a linear combination of the three 
most critical Euler buckling modes' shape {'^fc, k — 1,... ,3}, so that it reads as follows: 



where ^Si^^,k — 1, ... ,3j are three independent Gaussian random variates with mean 
jji^ — and standard deviation cr^ f=« 9.5 mm. This empirical model was determined with 
the following objectives: 

• the shape imperfection is Gaussian with a zero mean, so that the random variables 
should also be Gaussian with a zero mean; 

• the maximum amplitude of the shape imperfection around the mean surface (/>t^ = 0) 
±2(7^ is set equal to the half of the mean thickness jij^ — 76 mm, leading to cr^ ^ 



The other three random fields are assumed independent and lognormal with constant 
means pi^ = 76 mm, jJL^ — 200, 000 MPa and jji^^ — 390 MPa and coefficients of variation 

= 5%, = 3% and 5^^ — 7%. They are represented by three Karhunen-Loeve ex- 
pansions of three independent standard Gaussian random fields, whose sample paths are 
translated into lognormal sample paths - see Dubourg et al. (2009) for the mapping. These 
three Gaussian random fields are assumed to have the same isotropic squared exponential 
autocorrelation function with a correlation length i — 3, 500 mm. 

Due to the choice of this correlation function, the Fredholm integral equation involved in 
the Karhunen-Loeve discretization scheme has no analytical solution. A so-called wavelet- 
Galerkin numerical procedure was thus used as detailed in Phoon et al. (2002). Each 
random field is simulated by means of 30 independent standard Gaussian random variates 
leading to a relative mean squared discretization error of 3.70%. Finally the complete 
stochastic model involves 93 independent random variables. 

The reliability results are gathered in Table 2. The proposed importance sampling 
scheme leads to a failure probability in full agreement with the value obtained by sub- 
set simulation (as implemented in the FERUM toolbox v4.0 by Bourinet et al., 2009) which 



3 




(43) 



k=l 



9.5 mm. 
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Figure 6: One of the 4 most probable failure configurations. 



validates the proposed strategy. For this example the augmented failure probability pj ^ is 
equal to 2.06 x 10"'^ (with a coefficient of variation of 5.70%), and the correction factor 
Scoi-r is equal to 0.641 (with a coefficient of variation of 12.49%). 

A multiple FORM analysis revealed the existence of 4 most probable failure configu- 
rations identified by means of the restarted i-HLRF algorithm from Der Kiureghian and 
Dakessian (1998). These 4 failure modes corresponds to the 4 extreme cases for which the 
"demand" random field is maximal in the corner of the roof whereas the "capacity" ran- 
dom fields E, ay and h are minimal - one of these configurations is illustrated in Figure 6. 
The combination of these 4 failure modes in a serial system then allowed to give a Multi- 
FORM approximation of the failure probability (third column in Table 2) . In this case, the 
Ditlevsen's bounds (Ditlevsen and Madsen, 1996) coincide (up to the accuracy provided in 
Table 2). It thus proves the ability of the proposed strategy to deal with reliability problems 
featuring multiple design points in a reasonably high dimension. 

7. Conclusion 

Starting from the double premise that the usual surrogate-based reliability analyses do 
not permit to quantify the error made by using the metamodel instead of the original limit- 
state function, and that the existing variance reduction techniques remain time-consuming 
when the performance function involves the output of an expensive-to-evaluate black box 
function, an hybrid strategy has been proposed. 
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First, a probabilistic classification function based on the post-processing of a kriging pre- 
diction was introduced. This function allows a smoother classification than its deterministic 
counterpart (i.e. the indicator function of the failure domain) accounting for the epistemic 
uncertainty in the kriging prediction. The probabilistic classification is then used to for- 
mulate a quasi-optimal importance sampling density. Using elementary algebra the failure 
probability is recast as a product of two terms, namely the augmented failure probability Pf ^ 
which is evaluated by means of the meta-model only, and a correction factor a^orr that is 
computed from evaluations of the original limit-state function. In order to decide whether 
the kriging surrogate is accurate enough, a leave-one-out estimate of a^^^^ is used and the 
iterative refinement is stopped when it is in the order of magnitude of 1. Once the kriging 
surrogate has been built, the two terms of the product defining the failure probability may 
be evaluated in parallel. 

The method turned out to be efficient on several application examples, two of which 
are presented in this paper. It can handle problems featuring a reasonably high number of 
random variables and multiple design points. Further work is in progress to include the 
proposed algorithm within a reliability-based design optimization framework. 
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Appendix A. Calculation of the coefficient of variation 

The estimator in Eq. (32) is defined as the product of two estimators. Let us denote 
these two unbiased independent estimators as Pi, with variances cj^ for the sake of 
clarity. The calculation of the variance of the final estimator p — P1P2 proceeds as follows. 

First, according to its definition, the variance reads: 

al = Var [p] = E [p,^ P2'] - E [p.p^] ' (A.1) 
Since the two estimators p\ and p^ are independent, the variance also reads: 

o-H = E [p,^] E [p/] - E [p,] ' E [P2] ' (A.2) 
which may be further elaborated according to the Konig-Huyghens theorem: 

al = (e [p,] ' + al) (e [P2] ' + (jI) - E [p^] ' E [p^] ' (A.3) 
Due to the unbiasedness of the estimators, one finally gets: 

= {pI + ^D {pI + ^D-pIpI CA.4) 

= a^al+p^al+plal (A.5) 

Denoting by 5; = cjjpi the coefficients of variation of pi, i — 1, 2, the coefficient of 
variation oip — p\p2 eventually reads: 

5 = ^ = ^J 51 + 51 + 51 51 (A.6) 
P1P2 
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In practice usual target coefficients of variation 5(j,rget range from 1% to 10% so that 



5 \/5^ + 5l (A.7) 



Appendix B. The slice sampling technique 

This simulation technique consists in introducing an auxiliary scalar random variate U 
such that the joint PDF of the augmented vector {U ,X) reads: 

^. , / 1 if 0<u<p(a:) 

^^"'^^ = 1 otherwise ^^"^^ 

Note that this definition is nothing but a translation of the fundamental theorem of sim- 
ulation (see Robert and Casella, 2004, Theorem 2.15) whose key idea is to sample the 
augmented vector "under" the density curve (or surface) p(x) - see also Figure B.7. 

It is easy to show that the marginalization of J{u,x) w.r.t. u results in the target density 

P- 



J(u,x)du — 



ldu = p(x) (B.2) 



Hence, sampling X from p is equivalent to sampling (f/,X) from J and then ignoring U. 

The slice sampling procedure relies on the conditional distributions derived from J 
which read: 

U\X = x ~ {[0,p(x)]) (B.3) 
X|[/ = u ~ = {x eE" :p(x) > u}) (B.4) 

The slice y may not be easy to identify in high dimensions, but Neal (2003) proposed 
interesting algorithmic improvements which make the approach scalable to a wide variety 
of pseudo-PDFs. 

The two basic steps of the slice sampling algorithm are illustrated in Figure B.7. As 
other Markov chain sampling techniques, it does not require the explicit knowledge of 
the normalizing constant M of the target PDF p so that the value of M p{x) is sufficient. 
At iteration t, starting from a seed x^'^~^^ that is distributed according to the target PDF 
(i.e. such that M p{x'^'^~^^) / 0), one can first (i) generate a random number u^^^ which is 
uniformly distributed between and M p{x'^^^), and then (ii) generate a random vector x*-'-* 
that has uniform distribution on the t-th slice y'^^^ = jx e E" : Mp(x) > u"^'^]-. 
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